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

Last change on this file since 4339 was 4329, checked in by motisi, 5 years ago

Renamed wall_flags_0 to wall_flags_static_0

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