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

Last change on this file since 4472 was 4472, checked in by Giersch, 4 years ago

Profile output of the Kolmogorov length scale added

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