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

Last change on this file since 4571 was 4565, checked in by oliver.maas, 4 years ago

added new surface temperature forcing method for bc_pt_b = 'dirichlet': surface temperature pt(0) can be linearly increased by pt_surface_heating_rate (in K/h)

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