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

Last change on this file since 4670 was 4669, checked in by pavelkrc, 4 years ago

Fix multiple issues with radiation call times

  • Property svn:keywords set to Id
File size: 72.5 KB
Line 
1!> @file time_integration.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
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-2020 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: time_integration.f90 4669 2020-09-09 13:43:47Z maronga $
27! - Fix missing call of radiation after spinup
28! - Fix calculation of force_radiation_call
29! - Fix calculation of radiation times
30!
31! 4668 2020-09-09 13:00:16Z pavelkrc
32! Improve debug messages during timestepping
33!
34! 4581 2020-06-29 08:49:58Z suehring
35! Omit explicit pressure forcing via geostrophic wind components in case of
36! mesoscale nesting.
37!
38! 4578 2020-06-25 15:43:32Z gronemeier
39! bugfix: removed unused variables
40!
41! 4573 2020-06-24 13:08:47Z oliver.maas
42! calculate pt(0) incrementally by using dt_3d instead of calculating it absolutely
43! by using time_since_reference_point, because time_since_reference_point is set
44! to zero for initializing_actions = 'cyclic_fill'
45!
46! 4565 2020-06-15 08:30:38Z oliver.maas
47! added new surface temperature forcing method for bc_pt_b = 'dirichlet':
48! surface temperature pt(0) can be linearly increased by pt_surface_heating_rate (in K/h)
49!
50! 4564 2020-06-12 14:03:36Z raasch
51! Vertical nesting method of Huq et al. (2019) removed
52!
53! 4521 2020-05-06 11:39:49Z schwenkel
54! Rename variable
55!
56! 4511 2020-04-30 12:20:40Z raasch
57! chemistry decycling replaced by explicit setting of lateral boundary conditions
58!
59! 4508 2020-04-24 13:32:20Z raasch
60! salsa decycling replaced by explicit setting of lateral boundary conditions
61!
62! 4502 2020-04-17 16:14:16Z schwenkel
63! Implementation of ice microphysics
64!
65! 4472 2020-03-24 12:21:00Z Giersch
66! OPENACC COPYIN directive for ddx and ddy added
67!
68! 4466 2020-03-20 16:14:41Z suehring
69! Add advection fluxes to ACC copyin
70!
71! 4457 2020-03-11 14:20:43Z raasch
72! use statement for exchange horiz added
73!
74! 4444 2020-03-05 15:59:50Z raasch
75! bugfix: cpp-directives for serial mode added
76!
77! 4420 2020-02-24 14:13:56Z maronga
78! Added output control for wind turbine model
79!
80! 4403 2020-02-12 13:08:46Z banzhafs
81! Allowing both existing and on-demand emission read modes
82!
83! 4360 2020-01-07 11:25:50Z suehring
84! Bugfix, hour_call_emis uninitialized at first call of time_integration
85!
86! 4346 2019-12-18 11:55:56Z motisi
87! Introduction of wall_flags_total_0, which currently sets bits based on static
88! topography information used in wall_flags_static_0
89!
90! 4329 2019-12-10 15:46:36Z motisi
91! Renamed wall_flags_0 to wall_flags_static_0
92!
93! 4281 2019-10-29 15:15:39Z schwenkel
94! Moved boundary conditions to module interface
95!
96! 4276 2019-10-28 16:03:29Z schwenkel
97! Further modularization of lpm code components
98!
99! 4275 2019-10-28 15:34:55Z schwenkel
100! Move call oft lpm to the end of intermediate timestep loop
101!
102! 4268 2019-10-17 11:29:38Z schwenkel
103! Removing module specific boundary conditions an put them into their modules
104!
105! 4227 2019-09-10 18:04:34Z gronemeier
106! implement new palm_date_time_mod
107!
108! 4226 2019-09-10 17:03:24Z suehring
109! Changes in interface for the offline nesting
110!
111! 4182 2019-08-22 15:20:23Z scharf
112! Corrected "Former revisions" section
113!
114! 4170 2019-08-19 17:12:31Z gronemeier
115! copy diss, diss_p, tdiss_m to GPU
116!
117! 4144 2019-08-06 09:11:47Z raasch
118! relational operators .EQ., .NE., etc. replaced by ==, /=, etc.
119!
120! 4126 2019-07-30 11:09:11Z gronemeier
121! renamed routine to calculate uv exposure
122!
123! 4111 2019-07-22 18:16:57Z suehring
124! advc_flags_1 / advc_flags_2 renamed to advc_flags_m / advc_flags_s
125!
126! 4069 2019-07-01 14:05:51Z Giersch
127! Masked output running index mid has been introduced as a local variable to
128! avoid runtime error (Loop variable has been modified) in time_integration
129!
130! 4064 2019-07-01 05:33:33Z gronemeier
131! Moved call to radiation module out of intermediate time loop
132!
133! 4048 2019-06-21 21:00:21Z knoop
134! Moved production_e_init call into turbulence_closure_mod
135!
136! 4047 2019-06-21 18:58:09Z knoop
137! Added remainings of swap_timelevel upon its dissolution
138!
139! 4043 2019-06-18 16:59:00Z schwenkel
140! Further LPM modularization
141!
142! 4039 2019-06-18 10:32:41Z suehring
143! Rename subroutines in module for diagnostic quantities
144!
145! 4029 2019-06-14 14:04:35Z raasch
146! exchange of ghost points and boundary conditions separated for chemical species and SALSA module,
147! bugfix: decycling of chemistry species after nesting data transfer
148!
149! 4022 2019-06-12 11:52:39Z suehring
150! Call synthetic turbulence generator at last RK3 substep right after boundary
151! conditions are updated in offline nesting in order to assure that
152! perturbations are always imposed
153!
154! 4017 2019-06-06 12:16:46Z schwenkel
155! Mass (volume) flux correction included to ensure global mass conservation for child domains.
156!
157! 3994 2019-05-22 18:08:09Z suehring
158! output of turbulence intensity added
159!
160! 3988 2019-05-22 11:32:37Z kanani
161! Implement steerable output interval for virtual measurements
162!
163! 3968 2019-05-13 11:04:01Z suehring
164! replace nspec_out with n_matched_vars
165!
166! 3929 2019-04-24 12:52:08Z banzhafs
167! Reverse changes back from revision 3878: use chem_boundary_conds instead of
168! chem_boundary_conds_decycle
169!
170!
171! 3885 2019-04-11 11:29:34Z kanani
172! Changes related to global restructuring of location messages and introduction
173! of additional debug messages
174!
175! 3879 2019-04-08 20:25:23Z knoop
176! Moved wtm_forces to module_interface_actions
177!
178! 3872 2019-04-08 15:03:06Z knoop
179! Modifications made for salsa:
180! - Call salsa_emission_update at each time step but do the checks within
181!   salsa_emission_update (i.e. skip_time_do_salsa >= time_since_reference_point
182!   and next_aero_emission_update <= time_since_reference_point ).
183! - Renamed nbins --> nbins_aerosol, ncc_tot --> ncomponents_mass and
184!   ngast --> ngases_salsa and loop indices b, c and sg to ib, ic and ig
185! - Apply nesting for salsa variables
186! - Removed cpu_log calls speciffic for salsa.
187!
188! 3833 2019-03-28 15:04:04Z forkel
189! added USE chem_gasphase_mod, replaced nspec by nspec since fixed compounds are not integrated
190!
191! 3820 2019-03-27 11:53:41Z forkel
192! renamed do_emiss to emissions_anthropogenic (ecc)
193!
194!
195! 3774 2019-03-04 10:52:49Z moh.hefny
196! rephrase if statement to avoid unallocated array in case of
197! nesting_offline is false (crashing during debug mode)
198!
199! 3761 2019-02-25 15:31:42Z raasch $
200! module section re-formatted and openacc required variables moved to separate section,
201! re-formatting to 100 char line width
202!
203! 3745 2019-02-15 18:57:56Z suehring
204! Call indoor model after first timestep
205!
206! 3744 2019-02-15 18:38:58Z suehring
207! - Moved call of bio_calculate_thermal_index_maps from biometeorology module to
208! time_integration to make sure averaged input is updated before calculating.
209!
210! 3739 2019-02-13 08:05:17Z dom_dwd_user
211! Removed everything related to "time_bio_results" as this is never used.
212!
213! 3724 2019-02-06 16:28:23Z kanani
214! Correct double-used log_point_s unit
215!
216! 3719 2019-02-06 13:10:18Z kanani
217! - removed wind_turbine cpu measurement, since same time is measured inside
218!   wtm_forces subroutine as special measures
219! - moved the numerous vnest cpulog to special measures
220! - extended radiation cpulog over entire radiation part,
221!   moved radiation_interactions cpulog to special measures
222! - moved some cpu_log calls to this routine for better overview
223!
224! 3705 2019-01-29 19:56:39Z suehring
225! Data output for virtual measurements added
226!
227! 3704 2019-01-29 19:51:41Z suehring
228! Rename subroutines for surface-data output
229!
230! 3647 2019-01-02 14:10:44Z kanani
231! Bugfix: add time_since_reference_point to IF clause for data_output calls
232! (otherwise skip_time_* values don't come into affect with dt_do* = 0.0).
233! Clean up indoor_model and biometeorology model call.
234!
235! Revision 1.1  1997/08/11 06:19:04  raasch
236! Initial revision
237!
238!
239! Description:
240! ------------
241!> Integration in time of the model equations, statistical analysis and graphic
242!> output
243!------------------------------------------------------------------------------!
244 SUBROUTINE time_integration
245
246
247    USE advec_ws,                                                                                  &
248        ONLY:  ws_statistics
249
250    USE arrays_3d,                                                                                 &
251        ONLY:  diss, diss_p, dzu, e_p, nc_p, ni_p, nr_p, prho, pt, pt_p, pt_init, q, qc_p, qr_p,   &
252               q_init, q_p, qi_p, ref_state, rho_ocean, sa_p, s_p, tend, u, u_p, v, vpt, v_p, w_p
253
254#if defined( __parallel )  &&  ! defined( _OPENACC )
255    USE arrays_3d,                                                                                 &
256        ONLY:  e, nc, ni, nr, qc, qi, qr, s, w
257#endif
258
259    USE biometeorology_mod,                                                                        &
260        ONLY:  bio_calculate_thermal_index_maps, thermal_comfort, bio_calculate_uv_exposure,       &
261               uv_exposure
262
263    USE bulk_cloud_model_mod,                                                                      &
264        ONLY: bulk_cloud_model, calc_liquid_water_content, collision_turbulence,                   &
265              microphysics_ice_phase, microphysics_morrison, microphysics_seifert
266
267    USE calc_mean_profile_mod,                                                                     &
268        ONLY:  calc_mean_profile
269
270    USE chem_emissions_mod,                                                                        &
271        ONLY:  chem_emissions_setup, chem_emissions_update_on_demand
272
273    USE chem_gasphase_mod,                                                                         &
274        ONLY:  nvar
275
276    USE chem_modules,                                                                              &
277        ONLY:  bc_cs_t_val, chem_species, communicator_chem, emissions_anthropogenic,              &
278               emiss_read_legacy_mode, n_matched_vars
279
280    USE control_parameters,                                                                        &
281        ONLY:  advected_distance_x, advected_distance_y, air_chemistry, average_count_3d,          &
282               averaging_interval, averaging_interval_pr, bc_lr_cyc, bc_ns_cyc, bc_pt_t_val,       &
283               bc_q_t_val, biometeorology, call_psolver_at_all_substeps,  child_domain,            &
284               constant_flux_layer, constant_heatflux, create_disturbances,                        &
285               dopr_n, constant_diffusion, coupling_mode, coupling_start_time,                     &
286               current_timestep_number, debug_output_timestep, debug_string,                       &
287               disturbance_created, disturbance_energy_limit, dist_range,                          &
288               do_sum, dt_3d, dt_averaging_input, dt_averaging_input_pr, dt_coupling,              &
289               dt_data_output_av, dt_disturb, dt_do2d_xy, dt_do2d_xz, dt_do2d_yz, dt_do3d,         &
290               dt_domask,dt_dopts, dt_dopr, dt_dopr_listing, dt_dots, dt_run_control,              &
291               end_time, first_call_lpm, first_call_mas, galilei_transformation, humidity,         &
292               indoor_model, intermediate_timestep_count, intermediate_timestep_count_max,         &
293               land_surface, large_scale_forcing, loop_optimization, lsf_surf, lsf_vert, masks,    &
294               multi_agent_system_end, multi_agent_system_start, nesting_offline, neutral,         &
295               nr_timesteps_this_run, nudging, ocean_mode, passive_scalar, pt_reference,           &
296               pt_slope_offset, pt_surface_heating_rate,                                           &
297               random_heatflux, rans_tke_e, run_coupled, salsa,                                    &
298               simulated_time, simulated_time_chr, skip_time_do2d_xy, skip_time_do2d_xz,           &
299               skip_time_do2d_yz, skip_time_do3d, skip_time_domask, skip_time_dopr,                &
300               skip_time_data_output_av, sloping_surface, stop_dt, surface_output,                 &
301               terminate_coupled, terminate_run, timestep_scheme, time_coupling, time_do2d_xy,     &
302               time_do2d_xz, time_do2d_yz, time_do3d, time_domask, time_dopr, time_dopr_av,        &
303               time_dopr_listing, time_dopts, time_dosp, time_dosp_av, time_dots, time_do_av,      &
304               time_do_sla, time_disturb, time_run_control, time_since_reference_point,            &
305               timestep_count, turbulent_inflow, turbulent_outflow, urban_surface,                 &
306               use_initial_profile_as_reference, use_single_reference_value, u_gtrans, v_gtrans,   &
307               virtual_flight, virtual_measurement, ws_scheme_mom, ws_scheme_sca
308
309#if defined( __parallel )
310    USE control_parameters,                                                                        &
311        ONLY:  rans_mode
312#endif
313
314    USE cpulog,                                                                                    &
315        ONLY:  cpu_log, log_point, log_point_s
316
317    USE diagnostic_output_quantities_mod,                                                          &
318        ONLY:  doq_calculate,                                                                      &
319               timestep_number_at_prev_calc
320
321    USE exchange_horiz_mod,                                                                        &
322        ONLY:  exchange_horiz
323
324    USE flight_mod,                                                                                &
325        ONLY:  flight_measurement
326
327    USE grid_variables,                                                                            &
328        ONLY:  ddx, ddy
329
330    USE indices,                                                                                   &
331        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, nzb, nzt
332
333    USE indoor_model_mod,                                                                          &
334        ONLY:  dt_indoor, im_main_heatcool, time_indoor
335
336    USE interfaces
337
338    USE kinds
339
340    USE land_surface_model_mod,                                                                    &
341        ONLY:  lsm_boundary_condition, lsm_energy_balance, lsm_soil_model, skip_time_do_lsm
342
343    USE lagrangian_particle_model_mod,                                                             &
344        ONLY:  lpm_data_output_ptseries
345
346    USE lsf_nudging_mod,                                                                           &
347        ONLY:  calc_tnudge, ls_forcing_surf, ls_forcing_vert, nudge_ref
348
349    USE module_interface,                                                                          &
350        ONLY:  module_interface_actions, module_interface_swap_timelevel,                          &
351               module_interface_boundary_conditions
352
353    USE multi_agent_system_mod,                                                                    &
354        ONLY:  agents_active, multi_agent_system
355
356    USE nesting_offl_mod,                                                                          &
357        ONLY:  nesting_offl_bc,                                                                    &
358               nesting_offl_input,                                                                 &
359               nesting_offl_interpolation_factor,                                                  &
360               nesting_offl_mass_conservation
361
362    USE netcdf_data_input_mod,                                                                     &
363        ONLY:  chem_emis, chem_emis_att
364
365    USE ocean_mod,                                                                                 &
366        ONLY:  prho_reference
367
368    USE palm_date_time_mod,                                                                        &
369        ONLY:  get_date_time
370
371    USE particle_attributes,                                                                       &
372        ONLY:  particle_advection, particle_advection_start, use_sgs_for_particles, wang_kernel
373
374    USE pegrid
375
376#if defined( __parallel )
377    USE pmc_interface,                                                                             &
378        ONLY:  nested_run, nesting_mode, pmci_boundary_conds, pmci_datatrans, pmci_synchronize,    &
379        pmci_ensure_nest_mass_conservation, pmci_ensure_nest_mass_conservation_vertical,           &
380        pmci_set_swaplevel
381#endif
382
383    USE progress_bar,                                                                              &
384        ONLY:  finish_progress_bar, output_progress_bar
385
386    USE prognostic_equations_mod,                                                                  &
387        ONLY:  prognostic_equations_cache, prognostic_equations_vector
388
389    USE radiation_model_mod,                                                                       &
390        ONLY: dt_radiation, force_radiation_call, radiation, radiation_control,                    &
391              radiation_interaction, radiation_interactions, skip_time_do_radiation, time_radiation
392
393    USE salsa_mod,                                                                                 &
394        ONLY: aerosol_number, aerosol_mass, bc_am_t_val, bc_an_t_val, bc_gt_t_val,                 &
395              communicator_salsa, nbins_aerosol, ncomponents_mass, ngases_salsa,                   &
396              salsa_boundary_conditions, salsa_emission_update, salsa_gas, salsa_gases_from_chem,  &
397              skip_time_do_salsa
398
399    USE spectra_mod,                                                                               &
400        ONLY: average_count_sp, averaging_interval_sp, calc_spectra, dt_dosp, skip_time_dosp
401
402    USE statistics,                                                                                &
403        ONLY:  flow_statistics_called, hom, pr_palm, sums_ls_l
404
405
406    USE surface_layer_fluxes_mod,                                                                  &
407        ONLY:  surface_layer_fluxes
408
409    USE surface_data_output_mod,                                                                   &
410        ONLY:  average_count_surf, averaging_interval_surf, dt_dosurf, dt_dosurf_av,               &
411               surface_data_output, surface_data_output_averaging, skip_time_dosurf,               &
412               skip_time_dosurf_av, time_dosurf, time_dosurf_av
413
414    USE surface_mod,                                                                               &
415        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
416
417    USE synthetic_turbulence_generator_mod,                                                        &
418        ONLY:  dt_stg_call, dt_stg_adjust, parametrize_inflow_turbulence, stg_adjust, stg_main,    &
419               time_stg_adjust, time_stg_call, use_syn_turb_gen
420
421    USE turbulence_closure_mod,                                                                    &
422        ONLY:  tcm_diffusivities
423
424    USE urban_surface_mod,                                                                         &
425        ONLY:  usm_boundary_condition, usm_material_heat_model, usm_material_model,                &
426               usm_surface_energy_balance, usm_green_heat_model
427
428    USE virtual_measurement_mod,                                                                   &
429        ONLY:  dt_virtual_measurement,                                                             &
430               time_virtual_measurement,                                                           &
431               vm_data_output,                                                                     &
432               vm_sampling,                                                                        &
433               vm_time_start
434
435
436    USE wind_turbine_model_mod,                                                                    &
437        ONLY:  dt_data_output_wtm, time_wtm, wind_turbine, wtm_data_output
438
439#if defined( _OPENACC )
440    USE arrays_3d,                                                                                 &
441        ONLY:  d, dd2zu, ddzu, ddzw,                                                               &
442               diss_l_u,                                                                           &
443               diss_l_v,                                                                           &
444               diss_l_w,                                                                           &
445               diss_s_u,                                                                           &
446               diss_s_v,                                                                           &
447               diss_s_w,                                                                           &
448               drho_air, drho_air_zw, dzw, e,                                                      &
449               flux_l_u,                                                                           &
450               flux_l_v,                                                                           &
451               flux_l_w,                                                                           &
452               flux_s_u,                                                                           &
453               flux_s_v,                                                                           &
454               flux_s_w,                                                                           &
455               heatflux_output_conversion,                                                         &
456               kh, km, momentumflux_output_conversion, nc, ni, nr, p, ptdf_x, ptdf_y, qc, qi, qr, rdf,     &
457               rdf_sc, rho_air, rho_air_zw, s, tdiss_m, te_m, tpt_m, tu_m, tv_m, tw_m, ug, u_init, &
458               u_stokes_zu, vg, v_init, v_stokes_zu, w, zu
459
460    USE control_parameters,                                                                        &
461        ONLY:  tsc
462
463    USE indices,                                                                                   &
464        ONLY:  advc_flags_m, advc_flags_s, nyn, nyng, nys, nysg, nz, nzb_max, wall_flags_total_0
465
466    USE statistics,                                                                                &
467        ONLY:  rmask, statistic_regions, sums_l, sums_l_l, sums_us2_ws_l,                          &
468               sums_wsus_ws_l, sums_vs2_ws_l, sums_wsvs_ws_l, sums_ws2_ws_l, sums_wspts_ws_l,      &
469               sums_wsqs_ws_l, sums_wssas_ws_l, sums_wsqcs_ws_l, sums_wsqrs_ws_l, sums_wsncs_ws_l, &
470               sums_wsnrs_ws_l, sums_wsss_ws_l, weight_substep, sums_salsa_ws_l, sums_wsqis_ws_l,  &
471               sums_wsnis_ws_l
472
473    USE surface_mod,                                                                               &
474        ONLY:  bc_h, enter_surface_arrays, exit_surface_arrays
475#endif
476
477
478    IMPLICIT NONE
479
480    CHARACTER (LEN=9) ::  time_to_string   !<
481
482    INTEGER(iwp) ::  hour                !< hour of current time
483    INTEGER(iwp) ::  hour_call_emis = -1 !< last hour where emission was called
484    INTEGER(iwp) ::  ib                  !< index for aerosol size bins
485    INTEGER(iwp) ::  ic                  !< index for aerosol mass bins
486    INTEGER(iwp) ::  icc                 !< additional index for aerosol mass bins
487    INTEGER(iwp) ::  ig                  !< index for salsa gases
488    INTEGER(iwp) ::  mid                 !< masked output running index
489    INTEGER(iwp) ::  n                   !< loop counter for chemistry species
490
491    REAL(wp) ::  dt_3d_old  !< temporary storage of timestep to be used for
492                            !< steering of run control output interval
493    REAL(wp) ::  time_since_reference_point_save  !< original value of
494                                                  !< time_since_reference_point
495
496
497!
498!-- Copy data from arrays_3d
499!$ACC DATA &
500!$ACC COPY(d(nzb+1:nzt,nys:nyn,nxl:nxr)) &
501!$ACC COPY(diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
502!$ACC COPY(e(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
503!$ACC COPY(u(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
504!$ACC COPY(v(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
505!$ACC COPY(w(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
506!$ACC COPY(kh(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
507!$ACC COPY(km(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
508!$ACC COPY(p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
509!$ACC COPY(pt(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
510
511!$ACC DATA &
512!$ACC COPYIN(diss_l_u(0:nz+1,nys:nyn,0), flux_l_u(0:nz+1,nys:nyn,0)) &
513!$ACC COPYIN(diss_l_v(0:nz+1,nys:nyn,0), flux_l_v(0:nz+1,nys:nyn,0)) &
514!$ACC COPYIN(diss_l_w(0:nz+1,nys:nyn,0), flux_l_w(0:nz+1,nys:nyn,0)) &
515!$ACC COPYIN(diss_s_u(0:nz+1,0), flux_s_u(0:nz+1,0)) &
516!$ACC COPYIN(diss_s_v(0:nz+1,0), flux_s_v(0:nz+1,0)) &
517!$ACC COPYIN(diss_s_w(0:nz+1,0), flux_s_w(0:nz+1,0)) &
518!$ACC COPY(diss_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
519!$ACC COPY(e_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
520!$ACC COPY(u_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
521!$ACC COPY(v_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
522!$ACC COPY(w_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
523!$ACC COPY(pt_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
524!$ACC COPY(tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
525!$ACC COPY(tdiss_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
526!$ACC COPY(te_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
527!$ACC COPY(tu_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
528!$ACC COPY(tv_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
529!$ACC COPY(tw_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
530!$ACC COPY(tpt_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
531
532!$ACC DATA &
533!$ACC COPYIN(rho_air(nzb:nzt+1), drho_air(nzb:nzt+1)) &
534!$ACC COPYIN(rho_air_zw(nzb:nzt+1), drho_air_zw(nzb:nzt+1)) &
535!$ACC COPYIN(zu(nzb:nzt+1)) &
536!$ACC COPYIN(dzu(1:nzt+1), dzw(1:nzt+1)) &
537!$ACC COPYIN(ddzu(1:nzt+1), dd2zu(1:nzt)) &
538!$ACC COPYIN(ddzw(1:nzt+1)) &
539!$ACC COPYIN(heatflux_output_conversion(nzb:nzt+1)) &
540!$ACC COPYIN(momentumflux_output_conversion(nzb:nzt+1)) &
541!$ACC COPYIN(rdf(nzb+1:nzt), rdf_sc(nzb+1:nzt)) &
542!$ACC COPYIN(ptdf_x(nxlg:nxrg), ptdf_y(nysg:nyng)) &
543!$ACC COPYIN(ref_state(0:nz+1)) &
544!$ACC COPYIN(u_init(0:nz+1), v_init(0:nz+1)) &
545!$ACC COPYIN(u_stokes_zu(nzb:nzt+1), v_stokes_zu(nzb:nzt+1)) &
546!$ACC COPYIN(pt_init(0:nz+1)) &
547!$ACC COPYIN(ug(0:nz+1), vg(0:nz+1))
548
549!
550!-- Copy data from control_parameters
551!$ACC DATA &
552!$ACC COPYIN(tsc(1:5))
553
554!
555!-- Copy data from grid_variables
556!$ACC DATA &
557!$ACC COPYIN(ddx, ddy)
558
559!
560!-- Copy data from indices
561!$ACC DATA &
562!$ACC COPYIN(advc_flags_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
563!$ACC COPYIN(advc_flags_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
564!$ACC COPYIN(wall_flags_total_0(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
565
566!
567!-- Copy data from surface_mod
568!$ACC DATA &
569!$ACC COPYIN(bc_h(0:1)) &
570!$ACC COPYIN(bc_h(0)%i(1:bc_h(0)%ns)) &
571!$ACC COPYIN(bc_h(0)%j(1:bc_h(0)%ns)) &
572!$ACC COPYIN(bc_h(0)%k(1:bc_h(0)%ns)) &
573!$ACC COPYIN(bc_h(1)%i(1:bc_h(1)%ns)) &
574!$ACC COPYIN(bc_h(1)%j(1:bc_h(1)%ns)) &
575!$ACC COPYIN(bc_h(1)%k(1:bc_h(1)%ns))
576
577!
578!-- Copy data from statistics
579!$ACC DATA &
580!$ACC COPYIN(hom(0:nz+1,1:2,1:4,0)) &
581!$ACC COPYIN(rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions)) &
582!$ACC COPYIN(weight_substep(1:intermediate_timestep_count_max)) &
583!$ACC COPY(sums_l(nzb:nzt+1,1:pr_palm,0)) &
584!$ACC COPY(sums_l_l(nzb:nzt+1,0:statistic_regions,0)) &
585!$ACC COPY(sums_us2_ws_l(nzb:nzt+1,0)) &
586!$ACC COPY(sums_wsus_ws_l(nzb:nzt+1,0)) &
587!$ACC COPY(sums_vs2_ws_l(nzb:nzt+1,0)) &
588!$ACC COPY(sums_wsvs_ws_l(nzb:nzt+1,0)) &
589!$ACC COPY(sums_ws2_ws_l(nzb:nzt+1,0)) &
590!$ACC COPY(sums_wspts_ws_l(nzb:nzt+1,0)) &
591!$ACC COPY(sums_wssas_ws_l(nzb:nzt+1,0)) &
592!$ACC COPY(sums_wsqs_ws_l(nzb:nzt+1,0)) &
593!$ACC COPY(sums_wsqcs_ws_l(nzb:nzt+1,0)) &
594!$ACC COPY(sums_wsqis_ws_l(nzb:nzt+1,0)) &
595!$ACC COPY(sums_wsqrs_ws_l(nzb:nzt+1,0)) &
596!$ACC COPY(sums_wsncs_ws_l(nzb:nzt+1,0)) &
597!$ACC COPY(sums_wsnis_ws_l(nzb:nzt+1,0)) &
598!$ACC COPY(sums_wsnrs_ws_l(nzb:nzt+1,0)) &
599!$ACC COPY(sums_wsss_ws_l(nzb:nzt+1,0)) &
600!$ACC COPY(sums_salsa_ws_l(nzb:nzt+1,0))
601
602!
603!-- Next statement is to avoid compiler warnings about unused variables. Please
604!-- remove in case that you are using them. ddx and ddy need to be defined in
605!-- time_integration because of ACC COPYIN directive.
606    ddx = ddx
607    ddy = ddy
608
609#if defined( _OPENACC )
610    CALL enter_surface_arrays
611#endif
612
613!
614!-- At beginning determine the first time step
615    CALL timestep
616
617#if defined( __parallel )
618!
619!-- Synchronize the timestep in case of nested run.
620    IF ( nested_run )  THEN
621!
622!--    Synchronization by unifying the time step.
623!--    Global minimum of all time-steps is used for all.
624       CALL pmci_synchronize
625    ENDIF
626#endif
627
628!
629!-- Determine and print out the run control quantities before the first time
630!-- step of this run. For the initial run, some statistics (e.g. divergence)
631!-- need to be determined first --> CALL flow_statistics at the beginning of
632!-- run_control
633    CALL run_control
634!
635!-- Data exchange between coupled models in case that a call has been omitted
636!-- at the end of the previous run of a job chain.
637    IF ( coupling_mode /= 'uncoupled'  .AND.  run_coupled )  THEN
638!
639!--    In case of model termination initiated by the local model the coupler
640!--    must not be called because this would again cause an MPI hang.
641       DO WHILE ( time_coupling >= dt_coupling  .AND.  terminate_coupled == 0 )
642          CALL surface_coupler
643          time_coupling = time_coupling - dt_coupling
644       ENDDO
645       IF (time_coupling == 0.0_wp  .AND.  time_since_reference_point < dt_coupling )  THEN
646          time_coupling = time_since_reference_point
647       ENDIF
648    ENDIF
649
650    CALL location_message( 'atmosphere (and/or ocean) time-stepping', 'start' )
651
652!
653!-- Start of the time loop
654    DO  WHILE ( simulated_time < end_time  .AND.  .NOT. stop_dt  .AND. .NOT. terminate_run )
655
656       CALL cpu_log( log_point_s(10), 'timesteps', 'start' )
657
658       IF ( debug_output_timestep )  THEN
659           WRITE( debug_string, * ) 'time_integration', simulated_time
660           CALL debug_message( debug_string, 'start' )
661       ENDIF
662
663!
664!--    Determine ug, vg and w_subs in dependence on data from external file
665!--    LSF_DATA
666       IF ( large_scale_forcing .AND. lsf_vert )  THEN
667           CALL ls_forcing_vert ( simulated_time )
668           sums_ls_l = 0.0_wp
669       ENDIF
670
671!
672!--    Set pt_init and q_init to the current profiles taken from
673!--    NUDGING_DATA
674       IF ( nudging )  THEN
675           CALL nudge_ref ( simulated_time )
676!
677!--        Store temperature gradient at the top boundary for possible Neumann
678!--        boundary condition
679           bc_pt_t_val = ( pt_init(nzt+1) - pt_init(nzt) ) / dzu(nzt+1)
680           bc_q_t_val  = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1)
681           IF ( air_chemistry )  THEN
682              DO  n = 1, nvar
683                 bc_cs_t_val = (  chem_species(n)%conc_pr_init(nzt+1)                              &
684                                - chem_species(n)%conc_pr_init(nzt) )                              &
685                               / dzu(nzt+1)
686              ENDDO
687           ENDIF
688           IF ( salsa  .AND.  time_since_reference_point >= skip_time_do_salsa )  THEN
689              DO  ib = 1, nbins_aerosol
690                 bc_an_t_val = ( aerosol_number(ib)%init(nzt+1) - aerosol_number(ib)%init(nzt) ) / &
691                               dzu(nzt+1)
692                 DO  ic = 1, ncomponents_mass
693                    icc = ( ic - 1 ) * nbins_aerosol + ib
694                    bc_am_t_val = ( aerosol_mass(icc)%init(nzt+1) - aerosol_mass(icc)%init(nzt) ) /&
695                                  dzu(nzt+1)
696                 ENDDO
697              ENDDO
698              IF ( .NOT. salsa_gases_from_chem )  THEN
699                 DO  ig = 1, ngases_salsa
700                    bc_gt_t_val = ( salsa_gas(ig)%init(nzt+1) - salsa_gas(ig)%init(nzt) ) /        &
701                                  dzu(nzt+1)
702                 ENDDO
703              ENDIF
704           ENDIF
705       ENDIF
706!
707!--    Input of boundary data.
708       IF ( nesting_offline )  CALL nesting_offl_input
709!
710!--    Execute all other module actions routines
711       CALL module_interface_actions( 'before_timestep' )
712
713!
714!--    Start of intermediate step loop
715       intermediate_timestep_count = 0
716       DO  WHILE ( intermediate_timestep_count < intermediate_timestep_count_max )
717
718          intermediate_timestep_count = intermediate_timestep_count + 1
719
720!
721!--       Set the steering factors for the prognostic equations which depend
722!--       on the timestep scheme
723          CALL timestep_scheme_steering
724
725!
726!--       Calculate those variables needed in the tendency terms which need
727!--       global communication
728          IF ( .NOT. use_single_reference_value  .AND.  .NOT. use_initial_profile_as_reference )   &
729          THEN
730!
731!--          Horizontally averaged profiles to be used as reference state in
732!--          buoyancy terms (WARNING: only the respective last call of
733!--          calc_mean_profile defines the reference state!)
734             IF ( .NOT. neutral )  THEN
735                CALL calc_mean_profile( pt, 4 )
736                ref_state(:)  = hom(:,1,4,0) ! this is used in the buoyancy term
737             ENDIF
738             IF ( ocean_mode )  THEN
739                CALL calc_mean_profile( rho_ocean, 64 )
740                ref_state(:)  = hom(:,1,64,0)
741             ENDIF
742             IF ( humidity )  THEN
743                CALL calc_mean_profile( vpt, 44 )
744                ref_state(:)  = hom(:,1,44,0)
745             ENDIF
746!
747!--          Assure that ref_state does not become zero at any level
748!--          ( might be the case if a vertical level is completely occupied
749!--            with topography ).
750             ref_state = MERGE( MAXVAL(ref_state), ref_state, ref_state == 0.0_wp )
751          ENDIF
752
753          IF ( ( ws_scheme_mom .OR. ws_scheme_sca )  .AND.  intermediate_timestep_count == 1 )     &
754          THEN
755             CALL ws_statistics
756          ENDIF
757!
758!--       In case of nudging calculate current nudging time scale and horizontal
759!--       means of u, v, pt and q
760          IF ( nudging )  THEN
761             CALL calc_tnudge( simulated_time )
762             CALL calc_mean_profile( u, 1 )
763             CALL calc_mean_profile( v, 2 )
764             CALL calc_mean_profile( pt, 4 )
765             CALL calc_mean_profile( q, 41 )
766          ENDIF
767!
768!--       Execute all other module actions routunes
769          CALL module_interface_actions( 'before_prognostic_equations' )
770!
771!--       Solve the prognostic equations. A fast cache optimized version with
772!--       only one single loop is used in case of Piascek-Williams advection
773!--       scheme. NEC vector machines use a different version, because
774!--       in the other versions a good vectorization is prohibited due to
775!--       inlining problems.
776          IF ( loop_optimization == 'cache' )  THEN
777             CALL prognostic_equations_cache
778          ELSEIF ( loop_optimization == 'vector' )  THEN
779             CALL prognostic_equations_vector
780          ENDIF
781
782!
783!--       Movement of agents in multi agent system
784          IF ( agents_active  .AND.  time_since_reference_point >= multi_agent_system_start  .AND. &
785               time_since_reference_point <= multi_agent_system_end  .AND.                         &
786               intermediate_timestep_count == 1 )                                                  &
787          THEN
788             CALL multi_agent_system
789             first_call_mas = .FALSE.
790          ENDIF
791
792!
793!--       Exchange of ghost points (lateral boundary conditions)
794          CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'start' )
795
796          CALL exchange_horiz( u_p, nbgp )
797          CALL exchange_horiz( v_p, nbgp )
798          CALL exchange_horiz( w_p, nbgp )
799          CALL exchange_horiz( pt_p, nbgp )
800          IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e_p, nbgp )
801          IF ( rans_tke_e  .OR.  wang_kernel  .OR.  collision_turbulence                           &
802               .OR.  use_sgs_for_particles )  THEN
803             IF ( rans_tke_e )  THEN
804                CALL exchange_horiz( diss_p, nbgp )
805             ELSE
806                CALL exchange_horiz( diss, nbgp )
807             ENDIF
808          ENDIF
809          IF ( ocean_mode )  THEN
810             CALL exchange_horiz( sa_p, nbgp )
811             CALL exchange_horiz( rho_ocean, nbgp )
812             CALL exchange_horiz( prho, nbgp )
813          ENDIF
814          IF ( humidity )  THEN
815             CALL exchange_horiz( q_p, nbgp )
816             IF ( bulk_cloud_model .AND. microphysics_morrison )  THEN
817                CALL exchange_horiz( qc_p, nbgp )
818                CALL exchange_horiz( nc_p, nbgp )
819             ENDIF
820             IF ( bulk_cloud_model .AND. microphysics_seifert )  THEN
821                CALL exchange_horiz( qr_p, nbgp )
822                CALL exchange_horiz( nr_p, nbgp )
823             ENDIF
824             IF ( bulk_cloud_model .AND. microphysics_ice_phase )  THEN
825                CALL exchange_horiz( qi_p, nbgp )
826                CALL exchange_horiz( ni_p, nbgp )
827             ENDIF
828          ENDIF
829          IF ( passive_scalar )  CALL exchange_horiz( s_p, nbgp )
830          IF ( air_chemistry )  THEN
831             DO  n = 1, nvar
832                CALL exchange_horiz( chem_species(n)%conc_p, nbgp,                                 &
833                                     alternative_communicator = communicator_chem )
834             ENDDO
835          ENDIF
836
837          IF ( salsa  .AND.  time_since_reference_point >= skip_time_do_salsa )  THEN
838             DO  ib = 1, nbins_aerosol
839                CALL exchange_horiz( aerosol_number(ib)%conc_p, nbgp,                              &
840                                     alternative_communicator = communicator_salsa )
841                DO  ic = 1, ncomponents_mass
842                   icc = ( ic - 1 ) * nbins_aerosol + ib
843                   CALL exchange_horiz( aerosol_mass(icc)%conc_p, nbgp,                            &
844                                        alternative_communicator = communicator_salsa )
845                ENDDO
846             ENDDO
847             IF ( .NOT. salsa_gases_from_chem )  THEN
848                DO  ig = 1, ngases_salsa
849                   CALL exchange_horiz( salsa_gas(ig)%conc_p, nbgp,                                &
850                                        alternative_communicator = communicator_salsa )
851                ENDDO
852             ENDIF
853          ENDIF
854
855          CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'stop' )
856
857!
858!--       Boundary conditions for the prognostic quantities (except of the
859!--       velocities at the outflow in case of a non-cyclic lateral wall) and
860!--       boundary conditions for module-specific variables
861          CALL module_interface_boundary_conditions
862
863!
864!--       Incrementing timestep counter
865          timestep_count = timestep_count + 1
866
867          CALL cpu_log( log_point(28), 'swap_timelevel', 'start' )
868!
869!--       Set the swap level for all modules
870          CALL module_interface_swap_timelevel( MOD( timestep_count, 2) )
871
872#if defined( __parallel )
873!
874!--       Set the swap level for steering the pmc data transfer
875          IF ( nested_run )  CALL pmci_set_swaplevel( MOD( timestep_count, 2) + 1 )  !> @todo: why the +1 ?
876#endif
877
878          CALL cpu_log( log_point(28), 'swap_timelevel', 'stop' )
879
880#if defined( __parallel )
881          IF ( nested_run )  THEN
882
883             CALL cpu_log( log_point(60), 'nesting', 'start' )
884!
885!--          Domain nesting. The data transfer subroutines pmci_parent_datatrans
886!--          and pmci_child_datatrans are called inside the wrapper
887!--          subroutine pmci_datatrans according to the control parameters
888!--          nesting_mode and nesting_datatransfer_mode.
889!--          TO_DO: why is nesting_mode given as a parameter here?
890             CALL pmci_datatrans( nesting_mode )
891
892             IF ( TRIM( nesting_mode ) == 'two-way' .OR.  nesting_mode == 'vertical' )  THEN
893
894                CALL cpu_log( log_point_s(92), 'exchange-horiz-nest', 'start' )
895!
896!--             Exchange_horiz is needed for all parent-domains after the
897!--             anterpolation
898                CALL exchange_horiz( u, nbgp )
899                CALL exchange_horiz( v, nbgp )
900                CALL exchange_horiz( w, nbgp )
901                IF ( .NOT. neutral )  CALL exchange_horiz( pt, nbgp )
902
903                IF ( humidity )  THEN
904
905                   CALL exchange_horiz( q, nbgp )
906
907                   IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
908                       CALL exchange_horiz( qc, nbgp )
909                       CALL exchange_horiz( nc, nbgp )
910                   ENDIF
911                   IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
912                       CALL exchange_horiz( qr, nbgp )
913                       CALL exchange_horiz( nr, nbgp )
914                   ENDIF
915                   IF ( bulk_cloud_model  .AND.  microphysics_ice_phase )  THEN
916                      CALL exchange_horiz( qi, nbgp )
917                      CALL exchange_horiz( ni, nbgp )
918                   ENDIF
919
920                ENDIF
921
922                IF ( passive_scalar )  CALL exchange_horiz( s, nbgp )
923
924                IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e, nbgp )
925
926                IF ( .NOT. constant_diffusion  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
927                   CALL exchange_horiz( diss, nbgp )
928                ENDIF
929
930                IF ( air_chemistry )  THEN
931                   DO  n = 1, nvar
932                      CALL exchange_horiz( chem_species(n)%conc, nbgp,                             &
933                                           alternative_communicator = communicator_chem )
934                   ENDDO
935                ENDIF
936
937                IF ( salsa  .AND. time_since_reference_point >= skip_time_do_salsa )  THEN
938                   DO  ib = 1, nbins_aerosol
939                      CALL exchange_horiz( aerosol_number(ib)%conc, nbgp,                          &
940                                           alternative_communicator = communicator_salsa )
941                      DO  ic = 1, ncomponents_mass
942                         icc = ( ic - 1 ) * nbins_aerosol + ib
943                         CALL exchange_horiz( aerosol_mass(icc)%conc, nbgp,                        &
944                                              alternative_communicator = communicator_salsa )
945                      ENDDO
946                   ENDDO
947                   IF ( .NOT. salsa_gases_from_chem )  THEN
948                      DO  ig = 1, ngases_salsa
949                         CALL exchange_horiz( salsa_gas(ig)%conc, nbgp,                            &
950                                              alternative_communicator = communicator_salsa )
951                      ENDDO
952                   ENDIF
953                ENDIF
954                CALL cpu_log( log_point_s(92), 'exchange-horiz-nest', 'stop' )
955
956             ENDIF
957
958!
959!--          Set boundary conditions again after interpolation and anterpolation.
960             CALL pmci_boundary_conds
961
962             CALL cpu_log( log_point(60), 'nesting', 'stop' )
963
964          ENDIF
965#endif
966
967!
968!--       Temperature offset must be imposed at cyclic boundaries in x-direction
969!--       when a sloping surface is used
970          IF ( sloping_surface )  THEN
971             IF ( nxl ==  0 )  pt(:,:,nxlg:nxl-1) = pt(:,:,nxlg:nxl-1) - pt_slope_offset
972             IF ( nxr == nx )  pt(:,:,nxr+1:nxrg) = pt(:,:,nxr+1:nxrg) + pt_slope_offset
973          ENDIF
974
975!
976!--       Increase temperature pt(0) according to pt_surface_heating_rate (convert from K/h to K/s)
977          IF ( pt_surface_heating_rate /= 0.0_wp .AND. intermediate_timestep_count == 1 ) THEN
978             pt(0,:,:) = pt(0,:,:) + dt_3d * pt_surface_heating_rate / 3600.0_wp
979          ENDIF
980
981!
982!--       Impose a turbulent inflow using the recycling method
983          IF ( turbulent_inflow )  CALL inflow_turbulence
984
985!
986!--       Set values at outflow boundary using the special outflow condition
987          IF ( turbulent_outflow )  CALL outflow_turbulence
988
989!
990!--       Impose a random perturbation on the horizontal velocity field
991          IF ( create_disturbances  .AND.  ( call_psolver_at_all_substeps  .AND.                   &
992               intermediate_timestep_count == intermediate_timestep_count_max )                    &
993               .OR. ( .NOT. call_psolver_at_all_substeps  .AND.  intermediate_timestep_count == 1 ) ) &
994          THEN
995             time_disturb = time_disturb + dt_3d
996             IF ( time_disturb >= dt_disturb )  THEN
997                IF ( disturbance_energy_limit /= 0.0_wp  .AND.                                     &
998                     hom(nzb+5,1,pr_palm,0) < disturbance_energy_limit )  THEN
999                   CALL disturb_field( 'u', tend, u )
1000                   CALL disturb_field( 'v', tend, v )
1001                ELSEIF ( ( .NOT. bc_lr_cyc  .OR.  .NOT. bc_ns_cyc )                                &
1002                         .AND. .NOT. child_domain  .AND.  .NOT.  nesting_offline )                 &
1003                THEN
1004!
1005!--                Runs with a non-cyclic lateral wall need perturbations
1006!--                near the inflow throughout the whole simulation
1007                   dist_range = 1
1008                   CALL disturb_field( 'u', tend, u )
1009                   CALL disturb_field( 'v', tend, v )
1010                   dist_range = 0
1011                ENDIF
1012                time_disturb = time_disturb - dt_disturb
1013             ENDIF
1014          ENDIF
1015
1016!
1017!--       Map forcing data derived from larger scale model onto domain
1018!--       boundaries. Further, update geostrophic wind components.
1019          IF ( nesting_offline  .AND.  intermediate_timestep_count ==                              &
1020                                       intermediate_timestep_count_max  )  THEN
1021!--          Determine interpolation factor before boundary conditions and geostrophic wind
1022!--          is updated.
1023             CALL nesting_offl_interpolation_factor
1024             CALL nesting_offl_bc
1025!              CALL nesting_offl_geostrophic_wind
1026          ENDIF
1027!
1028!--       Impose a turbulent inflow using synthetic generated turbulence.
1029          IF ( use_syn_turb_gen  .AND.                                                             &
1030               intermediate_timestep_count == intermediate_timestep_count_max )  THEN
1031             CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
1032             CALL stg_main
1033             CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
1034          ENDIF
1035!
1036!--       Ensure mass conservation. This need to be done after imposing
1037!--       synthetic turbulence and top boundary condition for pressure is set to
1038!--       Neumann conditions.
1039!--       Is this also required in case of Dirichlet?
1040          IF ( nesting_offline )  CALL nesting_offl_mass_conservation
1041!
1042!--       Reduce the velocity divergence via the equation for perturbation
1043!--       pressure.
1044          IF ( intermediate_timestep_count == 1  .OR.  call_psolver_at_all_substeps )  THEN
1045
1046#if defined( __parallel )
1047!
1048!--          Mass (volume) flux correction to ensure global mass conservation for child domains.
1049             IF ( child_domain )  THEN
1050                IF ( nesting_mode == 'vertical' )  THEN
1051                   CALL pmci_ensure_nest_mass_conservation_vertical
1052                ELSE
1053                   CALL pmci_ensure_nest_mass_conservation
1054                ENDIF
1055             ENDIF
1056#endif
1057             CALL pres
1058
1059          ENDIF
1060!
1061!--       Particle transport/physics with the Lagrangian particle model
1062!--       (only once during intermediate steps, because it uses an Euler-step)
1063!--       ### particle model should be moved before prognostic_equations, in order
1064!--       to regard droplet interactions directly
1065
1066          CALL module_interface_actions( 'after_pressure_solver' )
1067!
1068!--       Interaction of droplets with temperature and mixing ratio.
1069!--       Droplet condensation and evaporation is calculated within
1070!--       advec_particles.
1071!
1072!--       If required, compute liquid water content
1073          IF ( bulk_cloud_model )  THEN
1074             CALL calc_liquid_water_content
1075          ENDIF
1076!
1077!--       If required, compute virtual potential temperature
1078          IF ( humidity )  THEN
1079             CALL compute_vpt
1080          ENDIF
1081
1082!
1083!--       Compute the diffusion quantities
1084          IF ( .NOT. constant_diffusion )  THEN
1085
1086!
1087!--          Determine surface fluxes shf and qsws and surface values
1088!--          pt_surface and q_surface in dependence on data from external
1089!--          file LSF_DATA respectively
1090             IF ( ( large_scale_forcing .AND. lsf_surf ) .AND.                                     &
1091                 intermediate_timestep_count == intermediate_timestep_count_max )                  &
1092             THEN
1093                CALL ls_forcing_surf( simulated_time )
1094             ENDIF
1095
1096!
1097!--          First the vertical (and horizontal) fluxes in the surface
1098!--          (constant flux) layer are computed
1099             IF ( constant_flux_layer )  THEN
1100                CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'start' )
1101                CALL surface_layer_fluxes
1102                CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'stop' )
1103             ENDIF
1104!
1105!--          If required, solve the energy balance for the surface and run soil
1106!--          model. Call for horizontal as well as vertical surfaces
1107             IF ( land_surface .AND. time_since_reference_point >= skip_time_do_lsm)  THEN
1108
1109                CALL cpu_log( log_point(54), 'land_surface', 'start' )
1110!
1111!--             Call for horizontal upward-facing surfaces
1112                CALL lsm_energy_balance( .TRUE., -1 )
1113                CALL lsm_soil_model( .TRUE., -1, .TRUE. )
1114!
1115!--             Call for northward-facing surfaces
1116                CALL lsm_energy_balance( .FALSE., 0 )
1117                CALL lsm_soil_model( .FALSE., 0, .TRUE. )
1118!
1119!--             Call for southward-facing surfaces
1120                CALL lsm_energy_balance( .FALSE., 1 )
1121                CALL lsm_soil_model( .FALSE., 1, .TRUE. )
1122!
1123!--             Call for eastward-facing surfaces
1124                CALL lsm_energy_balance( .FALSE., 2 )
1125                CALL lsm_soil_model( .FALSE., 2, .TRUE. )
1126!
1127!--             Call for westward-facing surfaces
1128                CALL lsm_energy_balance( .FALSE., 3 )
1129                CALL lsm_soil_model( .FALSE., 3, .TRUE. )
1130
1131!
1132!--             At the end, set boundary conditons for potential temperature
1133!--             and humidity after running the land-surface model. This
1134!--             might be important for the nesting, where arrays are transfered.
1135                CALL lsm_boundary_condition
1136
1137
1138                CALL cpu_log( log_point(54), 'land_surface', 'stop' )
1139             ENDIF
1140!
1141!--          If required, solve the energy balance for urban surfaces and run
1142!--          the material heat model
1143             IF (urban_surface) THEN
1144                CALL cpu_log( log_point(74), 'urban_surface', 'start' )
1145
1146                CALL usm_surface_energy_balance( .FALSE. )
1147                IF ( usm_material_model )  THEN
1148                   CALL usm_green_heat_model
1149                   CALL usm_material_heat_model ( .FALSE. )
1150                ENDIF
1151
1152!
1153!--             At the end, set boundary conditons for potential temperature
1154!--             and humidity after running the urban-surface model. This
1155!--             might be important for the nesting, where arrays are transfered.
1156                CALL usm_boundary_condition
1157
1158                CALL cpu_log( log_point(74), 'urban_surface', 'stop' )
1159             ENDIF
1160!
1161!--          Compute the diffusion coefficients
1162             CALL cpu_log( log_point(17), 'diffusivities', 'start' )
1163             IF ( .NOT. humidity ) THEN
1164                IF ( ocean_mode )  THEN
1165                   CALL tcm_diffusivities( prho, prho_reference )
1166                ELSE
1167                   CALL tcm_diffusivities( pt, pt_reference )
1168                ENDIF
1169             ELSE
1170                CALL tcm_diffusivities( vpt, pt_reference )
1171             ENDIF
1172             CALL cpu_log( log_point(17), 'diffusivities', 'stop' )
1173
1174          ENDIF
1175
1176       ENDDO   ! Intermediate step loop
1177
1178!
1179!--    Will be used at some point by flow_statistics.
1180       !$ACC UPDATE &
1181       !$ACC HOST(sums_l_l(nzb:nzt+1,0:statistic_regions,0)) &
1182       !$ACC HOST(sums_us2_ws_l(nzb:nzt+1,0)) &
1183       !$ACC HOST(sums_wsus_ws_l(nzb:nzt+1,0)) &
1184       !$ACC HOST(sums_vs2_ws_l(nzb:nzt+1,0)) &
1185       !$ACC HOST(sums_wsvs_ws_l(nzb:nzt+1,0)) &
1186       !$ACC HOST(sums_ws2_ws_l(nzb:nzt+1,0)) &
1187       !$ACC HOST(sums_wspts_ws_l(nzb:nzt+1,0)) &
1188       !$ACC HOST(sums_wssas_ws_l(nzb:nzt+1,0)) &
1189       !$ACC HOST(sums_wsqs_ws_l(nzb:nzt+1,0)) &
1190       !$ACC HOST(sums_wsqcs_ws_l(nzb:nzt+1,0)) &
1191       !$ACC HOST(sums_wsqis_ws_l(nzb:nzt+1,0)) &
1192       !$ACC HOST(sums_wsqrs_ws_l(nzb:nzt+1,0)) &
1193       !$ACC HOST(sums_wsncs_ws_l(nzb:nzt+1,0)) &
1194       !$ACC HOST(sums_wsnis_ws_l(nzb:nzt+1,0)) &
1195       !$ACC HOST(sums_wsnrs_ws_l(nzb:nzt+1,0)) &
1196       !$ACC HOST(sums_wsss_ws_l(nzb:nzt+1,0)) &
1197       !$ACC HOST(sums_salsa_ws_l(nzb:nzt+1,0))
1198
1199!
1200!--    If required, calculate radiative fluxes and heating rates
1201       IF ( radiation  .AND.  time_since_reference_point >= skip_time_do_radiation )  THEN
1202
1203          time_radiation = time_radiation + dt_3d
1204
1205          IF ( time_radiation >= dt_radiation  .OR.  force_radiation_call )  THEN
1206
1207             CALL cpu_log( log_point(50), 'radiation', 'start' )
1208
1209             IF ( time_radiation >= dt_radiation )  THEN
1210                time_radiation = time_radiation - dt_radiation
1211             ENDIF
1212
1213!
1214!--          Adjust the current time to the time step of the radiation model.
1215!--          Needed since radiation is pre-calculated and stored only on apparent
1216!--          solar positions
1217             time_since_reference_point_save = time_since_reference_point
1218             time_since_reference_point = time_since_reference_point -                            &
1219                                          MODULO(time_since_reference_point, dt_radiation)
1220
1221             CALL radiation_control
1222
1223             IF ( ( urban_surface  .OR.  land_surface )  .AND.  radiation_interactions )  THEN
1224                CALL cpu_log( log_point_s(46), 'radiation_interaction', 'start' )
1225                CALL radiation_interaction
1226                CALL cpu_log( log_point_s(46), 'radiation_interaction', 'stop' )
1227             ENDIF
1228
1229!
1230!--          Return the current time to its original value
1231             time_since_reference_point = time_since_reference_point_save
1232!
1233!--          Reset forcing of radiation call
1234             force_radiation_call = .FALSE.
1235
1236             CALL cpu_log( log_point(50), 'radiation', 'stop' )
1237
1238          ENDIF
1239       ENDIF
1240
1241
1242!
1243!-- 20200203 (ECC)
1244!-- allows for emission update mode in legacy mode as well as on-demand mode
1245!-- note that under on-demand mode emission update is no longer restricted to
1246!-- an hourly frequency, but whenever the simulation time corresponds to an
1247!-- inrement in emission timestamp value
1248
1249!
1250!-- If required, consider chemical emissions
1251
1252       IF  ( air_chemistry .AND. emissions_anthropogenic )  THEN
1253
1254          IF  ( emiss_read_legacy_mode )  THEN
1255!
1256!-- get hourly index and updates emission data when the hour is passed
1257
1258             CALL get_date_time( time_since_reference_point, hour=hour )
1259
1260             IF  ( hour_call_emis /= hour )   THEN
1261
1262                CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
1263                hour_call_emis = hour
1264
1265             ENDIF
1266
1267          ELSE
1268
1269             CALL chem_emissions_update_on_demand
1270
1271          ENDIF
1272
1273       ENDIF
1274
1275
1276!
1277!--    If required, consider aerosol emissions for the salsa model
1278       IF ( salsa )  THEN
1279!
1280!--       Call emission routine to update emissions if needed
1281          CALL salsa_emission_update
1282
1283       ENDIF
1284!
1285!--    If required, calculate indoor temperature, waste heat, heat flux
1286!--    through wall, etc.
1287!--    dt_indoor steers the frequency of the indoor model calculations.
1288!--    Note, at first timestep indoor model is called, in order to provide
1289!--    a waste heat flux.
1290       IF ( indoor_model )  THEN
1291
1292          time_indoor = time_indoor + dt_3d
1293
1294          IF ( time_indoor >= dt_indoor  .OR.  current_timestep_number == 0 )  THEN
1295
1296             time_indoor = time_indoor - dt_indoor
1297
1298             CALL cpu_log( log_point(76), 'indoor_model', 'start' )
1299             CALL im_main_heatcool
1300             CALL cpu_log( log_point(76), 'indoor_model', 'stop' )
1301
1302          ENDIF
1303       ENDIF
1304!
1305!--    Increase simulation time and output times
1306       nr_timesteps_this_run      = nr_timesteps_this_run + 1
1307       current_timestep_number    = current_timestep_number + 1
1308       simulated_time             = simulated_time   + dt_3d
1309       time_since_reference_point = simulated_time - coupling_start_time
1310       simulated_time_chr         = time_to_string( time_since_reference_point )
1311
1312
1313       IF ( time_since_reference_point >= skip_time_data_output_av )  THEN
1314          time_do_av         = time_do_av       + dt_3d
1315       ENDIF
1316       IF ( time_since_reference_point >= skip_time_do2d_xy )  THEN
1317          time_do2d_xy       = time_do2d_xy     + dt_3d
1318       ENDIF
1319       IF ( time_since_reference_point >= skip_time_do2d_xz )  THEN
1320          time_do2d_xz       = time_do2d_xz     + dt_3d
1321       ENDIF
1322       IF ( time_since_reference_point >= skip_time_do2d_yz )  THEN
1323          time_do2d_yz       = time_do2d_yz     + dt_3d
1324       ENDIF
1325       IF ( time_since_reference_point >= skip_time_do3d    )  THEN
1326          time_do3d          = time_do3d        + dt_3d
1327       ENDIF
1328       DO  mid = 1, masks
1329          IF ( time_since_reference_point >= skip_time_domask(mid) )  THEN
1330             time_domask(mid)= time_domask(mid) + dt_3d
1331          ENDIF
1332       ENDDO
1333       IF ( time_since_reference_point >= skip_time_dosp )  THEN
1334          time_dosp       = time_dosp        + dt_3d
1335       ENDIF
1336       time_dots          = time_dots        + dt_3d
1337       IF ( .NOT. first_call_lpm )  THEN
1338          time_dopts      = time_dopts       + dt_3d
1339       ENDIF
1340       IF ( time_since_reference_point >= skip_time_dopr )  THEN
1341          time_dopr       = time_dopr        + dt_3d
1342       ENDIF
1343       time_dopr_listing  = time_dopr_listing + dt_3d
1344       time_run_control   = time_run_control + dt_3d
1345!
1346!--    Increment time-counter for surface output
1347       IF ( surface_output )  THEN
1348          IF ( time_since_reference_point >= skip_time_dosurf )  THEN
1349             time_dosurf    = time_dosurf + dt_3d
1350          ENDIF
1351          IF ( time_since_reference_point >= skip_time_dosurf_av )  THEN
1352             time_dosurf_av = time_dosurf_av + dt_3d
1353          ENDIF
1354       ENDIF
1355!
1356!--    Increment time-counter for virtual measurements
1357       IF ( virtual_measurement  .AND.  vm_time_start <= time_since_reference_point )  THEN
1358          time_virtual_measurement = time_virtual_measurement + dt_3d
1359       ENDIF
1360
1361!
1362!--    Increment time-counter for wind turbine data output
1363       IF ( wind_turbine )  THEN
1364          time_wtm = time_wtm + dt_3d
1365       ENDIF
1366
1367!
1368!--    In case of synthetic turbulence generation and parametrized turbulence
1369!--    information, update the time counter and if required, adjust the
1370!--    STG to new atmospheric conditions.
1371       IF ( use_syn_turb_gen  )  THEN
1372          IF ( parametrize_inflow_turbulence )  THEN
1373             time_stg_adjust = time_stg_adjust + dt_3d
1374             IF ( time_stg_adjust >= dt_stg_adjust )  THEN
1375                CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
1376                CALL stg_adjust
1377                CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
1378             ENDIF
1379          ENDIF
1380          time_stg_call = time_stg_call + dt_3d
1381       ENDIF
1382
1383!
1384!--    Data exchange between coupled models
1385       IF ( coupling_mode /= 'uncoupled'  .AND.  run_coupled )  THEN
1386          time_coupling = time_coupling + dt_3d
1387
1388!
1389!--       In case of model termination initiated by the local model
1390!--       (terminate_coupled > 0), the coupler must be skipped because it would
1391!--       cause an MPI intercomminucation hang.
1392!--       If necessary, the coupler will be called at the beginning of the
1393!--       next restart run.
1394          DO WHILE ( time_coupling >= dt_coupling  .AND.  terminate_coupled == 0 )
1395             CALL surface_coupler
1396             time_coupling = time_coupling - dt_coupling
1397          ENDDO
1398       ENDIF
1399
1400!
1401!--    Biometeorology calculation of stationary thermal indices
1402!--    Todo (kanani): biometeorology needs own time_... treatment.
1403!--                   It might be that time_do2d_xy differs from time_do3d,
1404!--                   and then we might get trouble with the biomet output,
1405!--                   because we can have 2d and/or 3d biomet output!!
1406       IF ( biometeorology                                                                         &
1407            .AND. ( ( time_do3d >= dt_do3d  .AND.  time_since_reference_point >= skip_time_do3d )  &
1408                  .OR.                                                                             &
1409            ( time_do2d_xy >= dt_do2d_xy  .AND.  time_since_reference_point >= skip_time_do2d_xy ) &
1410                    ) )  THEN
1411!
1412!--       If required, do thermal comfort calculations
1413          IF ( thermal_comfort )  THEN
1414             CALL bio_calculate_thermal_index_maps ( .FALSE. )
1415          ENDIF
1416!
1417!--       If required, do UV exposure calculations
1418          IF ( uv_exposure )  THEN
1419             CALL bio_calculate_uv_exposure
1420          ENDIF
1421       ENDIF
1422
1423!
1424!--    Execute alle other module actions routunes
1425       CALL module_interface_actions( 'after_integration' )
1426
1427!
1428!--    If Galilei transformation is used, determine the distance that the
1429!--    model has moved so far
1430       IF ( galilei_transformation )  THEN
1431          advected_distance_x = advected_distance_x + u_gtrans * dt_3d
1432          advected_distance_y = advected_distance_y + v_gtrans * dt_3d
1433       ENDIF
1434
1435!
1436!--    Check, if restart is necessary (because cpu-time is expiring or
1437!--    because it is forced by user) and set stop flag
1438!--    This call is skipped if the remote model has already initiated a restart.
1439       IF ( .NOT. terminate_run )  CALL check_for_restart
1440
1441!
1442!--    Carry out statistical analysis and output at the requested output times.
1443!--    The MOD function is used for calculating the output time counters (like
1444!--    time_dopr) in order to regard a possible decrease of the output time
1445!--    interval in case of restart runs
1446
1447!
1448!--    Set a flag indicating that so far no statistics have been created
1449!--    for this time step
1450       flow_statistics_called = .FALSE.
1451
1452!
1453!--    If required, call flow_statistics for averaging in time
1454       IF ( averaging_interval_pr /= 0.0_wp  .AND.                                                 &
1455            ( dt_dopr - time_dopr ) <= averaging_interval_pr  .AND.                                &
1456            time_since_reference_point >= skip_time_dopr )  THEN
1457          time_dopr_av = time_dopr_av + dt_3d
1458          IF ( time_dopr_av >= dt_averaging_input_pr )  THEN
1459             do_sum = .TRUE.
1460             time_dopr_av = MOD( time_dopr_av, MAX( dt_averaging_input_pr, dt_3d ) )
1461          ENDIF
1462       ENDIF
1463       IF ( do_sum )  CALL flow_statistics
1464
1465!
1466!--    Sum-up 3d-arrays for later output of time-averaged 2d/3d/masked data
1467       IF ( averaging_interval /= 0.0_wp  .AND.                                                    &
1468            ( dt_data_output_av - time_do_av ) <= averaging_interval  .AND.                        &
1469            time_since_reference_point >= skip_time_data_output_av )                               &
1470       THEN
1471          time_do_sla = time_do_sla + dt_3d
1472          IF ( time_do_sla >= dt_averaging_input )  THEN
1473             IF ( current_timestep_number > timestep_number_at_prev_calc )                         &
1474                CALL doq_calculate
1475
1476             CALL sum_up_3d_data
1477             average_count_3d = average_count_3d + 1
1478             time_do_sla = MOD( time_do_sla, MAX( dt_averaging_input, dt_3d ) )
1479          ENDIF
1480       ENDIF
1481!
1482!--    Average surface data
1483       IF ( surface_output )  THEN
1484          IF ( averaging_interval_surf /= 0.0_wp                                                   &
1485                .AND.  ( dt_dosurf_av - time_dosurf_av ) <= averaging_interval_surf                &
1486                .AND.  time_since_reference_point >= skip_time_dosurf_av )  THEN
1487             IF ( time_dosurf_av >= dt_averaging_input )  THEN
1488                CALL surface_data_output_averaging
1489                average_count_surf = average_count_surf + 1
1490             ENDIF
1491          ENDIF
1492       ENDIF
1493
1494!
1495!--    Calculate spectra for time averaging
1496       IF ( averaging_interval_sp /= 0.0_wp  .AND. ( dt_dosp - time_dosp ) <= averaging_interval_sp&
1497            .AND.  time_since_reference_point >= skip_time_dosp )  THEN
1498          time_dosp_av = time_dosp_av + dt_3d
1499          IF ( time_dosp_av >= dt_averaging_input_pr )  THEN
1500             CALL calc_spectra
1501             time_dosp_av = MOD( time_dosp_av, MAX( dt_averaging_input_pr, dt_3d ) )
1502          ENDIF
1503       ENDIF
1504
1505!
1506!--    Call flight module and output data
1507       IF ( virtual_flight )  THEN
1508          CALL flight_measurement
1509          CALL data_output_flight
1510       ENDIF
1511!
1512!--    Take virtual measurements
1513       IF ( virtual_measurement  .AND.  time_virtual_measurement >= dt_virtual_measurement         &
1514                                 .AND.  vm_time_start <= time_since_reference_point )  THEN
1515          CALL vm_sampling
1516          CALL vm_data_output
1517          time_virtual_measurement = MOD(      time_virtual_measurement,                           &
1518                                          MAX( dt_virtual_measurement, dt_3d ) )
1519       ENDIF
1520
1521!
1522!--    Output wind turbine data
1523       IF ( wind_turbine  .AND.  time_wtm >= dt_data_output_wtm )  THEN
1524          CALL wtm_data_output
1525          time_wtm = MOD( time_wtm, MAX( dt_data_output_wtm, dt_3d ) )
1526       ENDIF
1527
1528!
1529!--    Profile output (ASCII) on file
1530       IF ( time_dopr_listing >= dt_dopr_listing )  THEN
1531          CALL print_1d
1532          time_dopr_listing = MOD( time_dopr_listing, MAX( dt_dopr_listing, dt_3d ) )
1533       ENDIF
1534
1535!
1536!--    Graphic output for PROFIL
1537       IF ( time_dopr >= dt_dopr  .AND.  time_since_reference_point >= skip_time_dopr )  THEN
1538          IF ( dopr_n /= 0 )  CALL data_output_profiles
1539          time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_3d ) )
1540          time_dopr_av = 0.0_wp    ! due to averaging (see above)
1541       ENDIF
1542
1543!
1544!--    Graphic output for time series
1545       IF ( time_dots >= dt_dots )  THEN
1546          CALL data_output_tseries
1547          time_dots = MOD( time_dots, MAX( dt_dots, dt_3d ) )
1548       ENDIF
1549
1550!
1551!--    Output of spectra (formatted for use with PROFIL), in case of no
1552!--    time averaging, spectra has to be calculated before
1553       IF ( time_dosp >= dt_dosp  .AND.  time_since_reference_point >= skip_time_dosp )  THEN
1554          IF ( average_count_sp == 0 )  CALL calc_spectra
1555          CALL data_output_spectra
1556          time_dosp = MOD( time_dosp, MAX( dt_dosp, dt_3d ) )
1557       ENDIF
1558
1559!
1560!--    2d-data output (cross-sections)
1561       IF ( time_do2d_xy >= dt_do2d_xy  .AND.  time_since_reference_point >= skip_time_do2d_xy )  THEN
1562          IF ( current_timestep_number > timestep_number_at_prev_calc )                            &
1563             CALL doq_calculate
1564
1565          CALL data_output_2d( 'xy', 0 )
1566          time_do2d_xy = MOD( time_do2d_xy, MAX( dt_do2d_xy, dt_3d ) )
1567       ENDIF
1568       IF ( time_do2d_xz >= dt_do2d_xz  .AND.  time_since_reference_point >= skip_time_do2d_xz )  THEN
1569          IF ( current_timestep_number > timestep_number_at_prev_calc )                            &
1570
1571             CALL doq_calculate
1572          CALL data_output_2d( 'xz', 0 )
1573          time_do2d_xz = MOD( time_do2d_xz, MAX( dt_do2d_xz, dt_3d ) )
1574       ENDIF
1575       IF ( time_do2d_yz >= dt_do2d_yz  .AND.  time_since_reference_point >= skip_time_do2d_yz )  THEN
1576          IF ( current_timestep_number > timestep_number_at_prev_calc )                            &
1577             CALL doq_calculate
1578
1579          CALL data_output_2d( 'yz', 0 )
1580          time_do2d_yz = MOD( time_do2d_yz, MAX( dt_do2d_yz, dt_3d ) )
1581       ENDIF
1582
1583!
1584!--    3d-data output (volume data)
1585       IF ( time_do3d >= dt_do3d  .AND.  time_since_reference_point >= skip_time_do3d )  THEN
1586          IF ( current_timestep_number > timestep_number_at_prev_calc )                            &
1587             CALL doq_calculate
1588
1589          CALL data_output_3d( 0 )
1590          time_do3d = MOD( time_do3d, MAX( dt_do3d, dt_3d ) )
1591       ENDIF
1592
1593!
1594!--    Masked data output
1595       DO  mid = 1, masks
1596          IF ( time_domask(mid) >= dt_domask(mid)                                                  &
1597               .AND.  time_since_reference_point >= skip_time_domask(mid) )  THEN
1598             IF ( current_timestep_number > timestep_number_at_prev_calc )                         &
1599                CALL doq_calculate
1600
1601             CALL data_output_mask( 0, mid )
1602             time_domask(mid) = MOD( time_domask(mid), MAX( dt_domask(mid), dt_3d ) )
1603          ENDIF
1604       ENDDO
1605
1606!
1607!--    Output of time-averaged 2d/3d/masked data
1608       IF ( time_do_av >= dt_data_output_av                                                        &
1609            .AND.  time_since_reference_point >= skip_time_data_output_av )  THEN
1610          CALL average_3d_data
1611!
1612!--       Udate thermal comfort indices based on updated averaged input
1613          IF ( biometeorology  .AND.  thermal_comfort )  THEN
1614             CALL bio_calculate_thermal_index_maps ( .TRUE. )
1615          ENDIF
1616          CALL data_output_2d( 'xy', 1 )
1617          CALL data_output_2d( 'xz', 1 )
1618          CALL data_output_2d( 'yz', 1 )
1619          CALL data_output_3d( 1 )
1620          DO  mid = 1, masks
1621             CALL data_output_mask( 1, mid )
1622          ENDDO
1623          time_do_av = MOD( time_do_av, MAX( dt_data_output_av, dt_3d ) )
1624       ENDIF
1625!
1626!--    Output of surface data, instantaneous and averaged data
1627       IF ( surface_output )  THEN
1628          IF ( time_dosurf >= dt_dosurf  .AND.  time_since_reference_point >= skip_time_dosurf )  THEN
1629             CALL surface_data_output( 0 )
1630             time_dosurf = MOD( time_dosurf, MAX( dt_dosurf, dt_3d ) )
1631          ENDIF
1632          IF ( time_dosurf_av >= dt_dosurf_av  .AND.  time_since_reference_point >= skip_time_dosurf_av )  THEN
1633             CALL surface_data_output( 1 )
1634             time_dosurf_av = MOD( time_dosurf_av, MAX( dt_dosurf_av, dt_3d ) )
1635          ENDIF
1636       ENDIF
1637
1638!
1639!--    Output of particle time series
1640       IF ( particle_advection )  THEN
1641          IF ( time_dopts >= dt_dopts  .OR.                                                        &
1642               ( time_since_reference_point >= particle_advection_start  .AND.                     &
1643                 first_call_lpm ) )  THEN
1644             CALL lpm_data_output_ptseries
1645             time_dopts = MOD( time_dopts, MAX( dt_dopts, dt_3d ) )
1646          ENDIF
1647       ENDIF
1648
1649!
1650!--    If required, set the heat flux for the next time step to a random value
1651       IF ( constant_heatflux  .AND.  random_heatflux )  THEN
1652          IF ( surf_def_h(0)%ns >= 1 )  THEN
1653             CALL cpu_log( log_point(23), 'disturb_heatflux', 'start' )
1654             CALL disturb_heatflux( surf_def_h(0) )
1655             CALL cpu_log( log_point(23), 'disturb_heatflux', 'stop' )
1656          ENDIF
1657          IF ( surf_lsm_h%ns    >= 1 )  THEN
1658             CALL cpu_log( log_point(23), 'disturb_heatflux', 'start' )
1659             CALL disturb_heatflux( surf_lsm_h    )
1660             CALL cpu_log( log_point(23), 'disturb_heatflux', 'stop' )
1661          ENDIF
1662          IF ( surf_usm_h%ns    >= 1 )  THEN
1663             CALL cpu_log( log_point(23), 'disturb_heatflux', 'start' )
1664             CALL disturb_heatflux( surf_usm_h    )
1665             CALL cpu_log( log_point(23), 'disturb_heatflux', 'stop' )
1666          ENDIF
1667       ENDIF
1668
1669!
1670!--    Execute alle other module actions routunes
1671       CALL module_interface_actions( 'after_timestep' )
1672
1673!
1674!--    Determine size of next time step. Save timestep dt_3d because it is
1675!--    newly calculated in routine timestep, but required further below for
1676!--    steering the run control output interval
1677       dt_3d_old = dt_3d
1678       CALL timestep
1679
1680#if defined( __parallel )
1681!
1682!--    Synchronize the timestep in case of nested run.
1683       IF ( nested_run )  THEN
1684!
1685!--       Synchronize by unifying the time step.
1686!--       Global minimum of all time-steps is used for all.
1687          CALL pmci_synchronize
1688       ENDIF
1689#endif
1690
1691!
1692!--    Computation and output of run control parameters.
1693!--    This is also done whenever perturbations have been imposed
1694       IF ( time_run_control >= dt_run_control  .OR.                                               &
1695            timestep_scheme(1:5) /= 'runge'  .OR.  disturbance_created )                           &
1696       THEN
1697          CALL run_control
1698          IF ( time_run_control >= dt_run_control )  THEN
1699             time_run_control = MOD( time_run_control, MAX( dt_run_control, dt_3d_old ) )
1700          ENDIF
1701       ENDIF
1702
1703!
1704!--    Output elapsed simulated time in form of a progress bar on stdout
1705       IF ( myid == 0 )  CALL output_progress_bar
1706
1707       IF ( debug_output_timestep )  THEN
1708           WRITE( debug_string, * ) 'time_integration', simulated_time
1709           CALL debug_message( debug_string, 'end' )
1710       ENDIF
1711
1712       CALL cpu_log( log_point_s(10), 'timesteps', 'stop' )
1713
1714
1715    ENDDO   ! time loop
1716
1717#if defined( _OPENACC )
1718    CALL exit_surface_arrays
1719#endif
1720!$ACC END DATA
1721!$ACC END DATA
1722!$ACC END DATA
1723!$ACC END DATA
1724!$ACC END DATA
1725!$ACC END DATA
1726!$ACC END DATA
1727!$ACC END DATA
1728
1729    IF ( myid == 0 )  CALL finish_progress_bar
1730
1731    CALL location_message( 'atmosphere (and/or ocean) time-stepping', 'finished' )
1732
1733 END SUBROUTINE time_integration
Note: See TracBrowser for help on using the repository browser.