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

Last change on this file since 4578 was 4578, checked in by gronemeier, 7 months ago

message.f90:

  • bugfix : do not save input values from last call of routines debug_message and location_message
  • changes: layout changes according to PALM coding standards

time_integration.f90:

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