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

Last change on this file since 4878 was 4848, checked in by gronemeier, 4 years ago

bugfix: removed syn_turb_gen from restart files, replaced use_syn_turb_gen by syn_turb_gen

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