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

Last change on this file since 4403 was 4403, checked in by banzhafs, 20 months ago

chemistry model: implemented on-demand emission read mode for LOD 2

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