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

Last change on this file since 4828 was 4828, checked in by Giersch, 4 years ago

Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code

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