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

Last change on this file since 4287 was 4281, checked in by schwenkel, 4 years ago

Moved boundary_conds to dynamics module

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