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

Last change on this file since 4420 was 4420, checked in by maronga, 5 years ago

added steering for NetCDF output for wind turbine model; minor fix in palmrungui

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