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

Last change on this file since 4227 was 4227, checked in by gronemeier, 5 years ago

implement new palm_date_time_mod; replaced namelist parameters time_utc_init and day_of_year_init by origin_date_time

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