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

Last change on this file since 4444 was 4444, checked in by raasch, 3 years ago

bugfix: cpp-directives for serial mode added

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