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

Last change on this file since 4226 was 4226, checked in by suehring, 2 years ago

Offline nesting: data input modularized; time variable is defined relative to time_utc_init, so that input data is correctly mapped to actual model time; checks rephrased and new checks for the time dimension added; Netcdf input: routine to retrieve dimension length renamed

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