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

Last change on this file since 4782 was 4732, checked in by schwenkel, 4 years ago

add use statements for openacc

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