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

Last change on this file since 4356 was 4356, checked in by suehring, 5 years ago

Bugfix in message calls for local checks; error messages in init_grid slightly revised; bugfix in time_integration (uninitialized emission time index); read_restart_data (change from J.Resler): use dynamic array allocation rather than automatic arrays to avoid problems with stack memory

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