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

Last change on this file since 4268 was 4268, checked in by schwenkel, 4 years ago

Introducing module interface for boundary conditions and move module specific boundary conditions into their modules

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