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

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