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

Last change on this file since 4346 was 4346, checked in by motisi, 16 months ago

Introduction of wall_flags_total_0, which currently sets bits based on static topography information used in wall_flags_static_0

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