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

Last change on this file since 4466 was 4466, checked in by suehring, 4 years ago

Vector branch in advec_ws optimized, symmetric boundary conditions implemented in vector branch

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