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

Last change on this file since 4581 was 4581, checked in by suehring, 12 months ago

mesoscale nesting: omit explicit pressure forcing via geostrophic wind components; chemistry: enable profile output of vertical fluxes; urban-surface: bugfix in initialization in case of cyclic_fill

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