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

Last change on this file since 4276 was 4276, checked in by schwenkel, 18 months ago

modularize lpm code components of time integration

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