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

Last change on this file since 4725 was 4671, checked in by pavelkrc, 4 years ago

Radiative transfer model RTM version 4.1

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